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Abstract 

Pressure solution is an important process in 
sedimentary basins, and its behaviour depends 
mainly on the sediment rheology and temperature 
distribution. The compaction relation of pressure 
solution is typically assumed to be a viscous one 
and is often written as a relationship between 
effective stress and strain rate. A new derivation 
of viscous compaction relation is formulated based 
on more realistic boundary conditions at grain 
contacts. A nonlinear diffusion problem with a 
moving boundary is solved numerically and a 
simple asymptotic solution is given to compare 
with numerical simulations. Pressure solution 
is significantly influenced by the temperature 
gradient. Porosity reduction due to pressure 
solution is enhanced in an environment with a 
higher thermal gradient, while porosity decreases 
much slowly in the region where the thermal 
gradient is small. Pressure solution tends to 
complete more quickly at shallower depths and 
earlier time in higher temperature environment 
than that in a low one. These features of pressure 
solution in porous sediments are analysed using 
a perturbation method to get a solution for the 
steady state. Comparison with real data shows a 
reasonably very good agreement. 
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1 INTRODUCTION 

Pressure solution is a very common and im- 
portant deformation process in porous media 



and granular materials such as sediments and 
soils. Pressure solution also occurs in sedimen- 
tary basins where hydrocarbons and oil are pri- 
marily formed. The modelling of such com- 
pactional flow is thus important in the oil in- 
dustry as well as in civil engineering. One par- 
ticular problem which affects drilling process 
is the occasional occurrence of abnormally high 
pore fluid pressures, which, if encountered sud- 
denly, can cause drill hole collapse and conse- 
quent failure of the drilling operation. There- 
fore, an industrially important objective is to 
predict over pressuring before drilling and to 
identify its precursors during drilling. An es- 
sential step to achieve such objectives is the sci- 
entific understanding of their mechanisms and 
the evolutionary history of post-depositional 
sediments such as shales. 

Compaction is the process of volume reduc- 
tion via pore- water expulsion within sediments 
due to the increasing weight of overburden 
load. The requirement of its occurrence is not 
only the application of an overburden load but 
also the expulsion of pore water. The extent of 
compaction is strongly influenced by sedimen- 
tation history and the lithology of sediments. 
The freshly deposited loosely packed sediments 
tend to evolve, like an open system, towards a 
closely packed grain framework during the ini- 
tial stages of burial compaction and this is ac- 
complished by the processes of grain slippage, 
rotation, bending and brittle fracturing. Such 
reorientation processes are collectively referred 
to as mechanical compaction, which generally 
takes place in the first 1-2 km of burial. Af- 
ter this initial porosity loss, further porosity 
reduction is accomplished by the process of 
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chemical compaction such as pressure solution 
at grain contacts [1,2,3]. 

Pressure solution has been considered as an 
important process in deformation and poros- 
ity change during compaction in sedimentary 
rocks [4,5]. Pressure solution refers to a pro- 
cess by which grains dissolve at intergranular 
contacts under non-hydrostatic stress and re- 
precipitate in pore spaces, thus resulting in 
compaction. The solubility of minerals in- 
creases with increasing effective stress at grain 
contacts. Pressure dissolution at grain con- 
tacts is therefore a compactional response of 
the sediment during burial in an attempt to in- 
crease the grain contact area so as to distribute 
the effective stress over a larger surface. Such 
a compaction process is typically assumed to 
viscous [5,6,7] and it is usually referred to as 
viscous compaction, viscous creep or pressure 
solution creep. Its rheological constitutive re- 
lation (or compaction relation) is often written 
as a relationship between effective stress and 
strain rate. 

A typical form of pressure solution is inter- 
granular pressure solution (IPS) which occurs 
at individual grain contacts and free face pres- 
sure solution (FFPS) which occurs at the face 
in contact with the pore fluid, but most stud- 
ies have concentrated on the former one (IPS). 
Extensive studies [1,5,6,7,8] on pressure solutin 
have been carried out in the last two decades, 
and a comprehensive literature review on these 
models was given by Tada and Siever [8]. A 
more recent and brief review can be found in 
[5,6]. Despite of its geological importance, the 
mechanism of pressure solution is still poorly 
understood. Recently, Fowler and Yang [5] 
present a new mathematical approach to model 
pressure solution and viscous compaction in 
sedimentary basins and show that the main pa- 
rameter controlling the compaction and poros- 
ity reduction is the compaction parameter A, 
the ratio of hydraulic conductivity to the sed- 
imentation rate. Compaction relation of poro- 
elastic and viscous type is also an important 
factor controlling the behaviour of compaction 
profile. However, the temperature effect has 
not been included in their approach. Thus, we 
mainly investigate the effect of different tem- 



perature gradients on the viscous compaction 
due to pressure solution in sedimentary basins. 

2 MATHEMATICAL 
MODEL 

For the convenience of investigating the effect 
of compaction in porous media due to pure 
density differences, we will assume the basic 
model of compaction is rather analogous to the 
process of soil consolidation. The porous me- 
dia act compressible porous matrix, so 
that mass conservation of pore fluid together 
with Darcy's law leads to the 1-D model equa- 
tions of the general type [1,5]. 
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(force balance), (4) 

where u l and u s are the velocities of fluid and 
solid matrix, k and p are the matrix perme- 
ability and the liquid viscosity, pi and p s are 
the densities of fluid and solid matrix, p is 
pore pressure, p e is the effective pressure, rj 
is medium viscosity and £ is compaction vis- 
cosity, and g is the gravitational acceleration. 
Combining the force balance and Darcy's law 
to eliminate p, we have 

(u i_ u .) = M { ^[(i + 0)p e ]_( p- _p l )(i_^ }> 

(5) 

which is a derived form of Darcy's law. By 
assuming the densities p s and pi are constants, 
we can see that only the density difference p s — 
pi is important to the flow evolution. Thus, the 
compactional flow is essentially density-driven 
flow in porous media. 

Compaction relation is a relationship be- 
tween effective pressure p e and strain rate e = 
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or porosity <fi [5,6,7]. The common ap- 
proach in soil mechanics and sediment com- 
paction is to model this generally nonlinear 
behaviour as poroelastic, that is to say, a rela- 
tionship of Athy's law type p e = p e (<f>), which 
is derived from fitting the real data of sedi- 
ments. However, this poroelastic compaction 
law is only valid for the compaction in porous 
media in the upper and shallow region, where 
compaction occurs due to the pure mechani- 
cal movements such as grain sliding and pack- 
ing rearrangement. In the more deeper region, 
mechanical compaction is gradually replaced 
by the chemical compaction due to stress- 
enhanced flow along the grain boundary from 
the grain contact areas to the free pore, where 
pressure is essentially pore pressure. A typical 
process of such chemical compaction in sedi- 
ment is pressure solution whose rheological be- 
havior is usually viscous, so that it sometimes 
called viscous pressure solution. 

The mathematical formulation for viscous 
compaction is to derive a relation between 
creep rate e and effective stress a e . Rutter's 
creep relation is widely used [7,8,9] 

_ A k c wD gb 
6 ~ p s d? ae ' W 

where a e is the effective normal stress across 
the grain contacts, A/, is a constant, cq is the 
equilibrium concentration (of quartz) in pore 
fluid, p, d are the density and (averaged) grain 
diameter (of quartz). D gb is the diffusivity of 
the solute in water along grain boundaries with 
a thickness w. D gb also varies with tempera- 
ture T 

D gb (T) = D gb e-^, (7) 

where E a is the effective activation energy with 
a value of 3 ~ 6 kJ/mole or even much lower [1, 
6] . From the values of the diffusion coefficient 
in quartz-water and rocksalt-water systems at 
300, 600, 1200 K, we get an estimate value of 
E a fa 0.65 kcal/mole [1,7]. 

Note that a e = -(1 + ||)p e and e = 
With this, (6) becomes the following com- 
paction law 

p e = -£V.u s . (8) 

More generally speaking, £ is also a function 
of porosity <f>. The compaction law is analo- 



gous to Fowler's viscous compaction laws used 
in studies of magma transport in the Earth's 
mantle. 

2.1 Derivation of Viscous Law 

The approach of deriving the law of viscous 
compaction depends on the underlying mech- 
anism. The classical theoretical considera- 
tion assumed a grain-boundary diffusion film of 
constant thickness and diffusivity, while others 
used the concept of a roughened, fluid-invaded 
non-equilibrium contact structure Shimizu [6] 
presented a kinetic approach extending Coble's 
classical treatment of grain boundary diffu- 
sion creep by including the kinetics of quartz 
dissolution/precipitation reaction. Shimuzu's 
derivation is instructive although the boundary 
conditions used in his formulation are question- 
able and unrealistic. In addition, Shimuzu's 
1-D approximation is only valid for a closed 
system due to §^ = used in his work when 
the thickness w of the water film is small with 
respect to the grain diameter (d) [5,10]. In or- 
der to correctly formulate the derivation, we 
now provide a new derivation by using more 
realistic boundary conditions in an open sys- 
tem. 

Now let us consider the intergranular contact 
region as a disk with a radius r = L. Let J(r) 
be the radial component of solute mass flux, e 
be the average strain rate, and v is the uniform 
shortening velocity of the upper grain relative 
to the lower grain due to the pressure solution 
creep [6,10]. The kinetic relation between v 
and e becomes 

v = ed. (9) 

For simplicity, we assume that the film thick- 
ness w is constant and the diffusion is near 
steady-state. Mass conservation gives 

2vrr J(r) + p s irr 2 v = 0, (10) 

where the flux J(r) obeys Fick's Law 

dc 

J(r) = ~D gb w-. (11) 

The steady-state solution of concentration c(r) 
for the boundary conditions c r = at r = 0, 



3 



c = Co at r = L is 

c(r) = c - 



(L 2 -r 2 ). (12) 



The parabolic change of concentration c(r) im- 
plies that the stress a(r) should be hetero- 
geneously distributed in the contact region. 
From a relation of effective stress and concen- 
tration [10], we have 



e, n R T , c(r) 

a (r) = ln^A 

v m c 



(13) 



where v m is the molar volume of the sediment. 
We have used here the condition a e (r) = at 
r = L. Let a be the averaged effective stress, 
then 



nL 2 a = / 2ira e (r)rdr. 
Jo 

Combining (13) and (14), we have 



(14) 



a 



2RT 

v m L 2 



f L rln[l 
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-{L 2 — r 2 )]dr. 



4c D gb w 

Using (9) and integrating by parts, we have 
1 
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where 



- — [(1- 
v m 11 BL 2 



B 



(15) 

we ha 1 
)ln(l - BL 2 ) - 1], (16) 

p s ed 



(17) 

ic D gb w 

By defining a critical effective stress a c (and 
equivalently a critical creep rate e c ) when 
BL 2 = 1 



Or. = 



RT . _ 4coD g bW 



e r . = 



p s L 2 d 

equation (16) can be rewritten as 



- = [l-(l--)ln(l--)]- 

a c e e c 



(18) 



(19) 



From the typical values of T ~ 300 K, R ~ 8.31 
J mol -1 K -1 , and v m ~ 2.6 x 10~ 5 m 3 mol -1 
[11], we can use the definition (18) to estimate 
the typical value of a Cl which is about 95 MPa. 
Clearly, if | a \ <^ a c , we have 

Av m coD gb w Wv m c D gb w 
e = „„ — -t^-cf = — — =k — c, (20) 



which is exactly the creep law. Here we have 
used L = d/2. A different choice of L = O(d) 
will only introduce an additional shape factor 
into the above relation. Under upper-crustal 
stress conditions a < 100 MPa, the above ap- 
proximation is valid as we expected. At higher 
stress states, we can use \a \ ^$> a c , then (19) 
becomes 



4c D gb w 
p s dL 2 



[1-e 



RT 



(21) 



Let L 2 = A(P/a s , and a s = O(l) is a shape 
factor. The above relation (21) becomes 



. _ a s c D gb w 
Psd 4 



(22) 



RTp s dL 2 



RTp s d 3 



which degenerates into (20) when v m ojRT 
1, but (20) may be inaccurate when | a \ ~ a c . 
The new compaction relation (22) is more ac- 
curate and valid in a more wide range of pa- 
rameter variations. 

Furthermore, the newly derived viscous 
compaction law (22) shows that the strain rate 
due to pressure solution is controlled by many 
parameters such as grain size (d) , grain geome- 
try (ctg), temperature (T), grain-boundary dif- 
fusion coeffient (D gb ). This is consistent with 
Dewers and Hajash's [12] empirical law de- 
rived from a quartz compaction experiment. 
However, since the complicated dependence 
on many parameters and nonlinear features in 
(22), various simplified version or approximate 
forms have been used by many authors in ear- 
lier work [13-16]. One common simplification 
of (22) is its linearised form such as (6) used by 
Rutter [7] and Ortoleva [13]. A slightly differ- 
ent formulation of this compaction law is ex- 
pressed in terms of porosity strain versus ef- 
fective stress (instead of using the strain rate). 
Schneider et al. [14] used a relationship be- 
tween porosity and effective stress while Lan- 
der and Walderhaug [15] used an exponential 
form of intergranular volume function of 
effective stress. Revil [16] used a relationship 
between porosity strain and effective stress, 
which includes time t explicitly in his formu- 
lation. However, these different formulations 
can be transformed into a relationship similar 
to (20) but such a transformation may depend 
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on the grain packing structure because of the 
calculation of porosity strain and porosity. For 
simplicity, we will only use the form (20) in the 
rest of the paper. 

2.2 Boundary conditions 

The boundary conditions for the governing 
equations follows. The bottom bound- 

ary at z = is assumed to be impermeable 

u s = u l = 0, (23) 

and a top condition at z = h is kinetic 

h = m s + u s , (24) 

where m s is the sedimentation rate at z = h. 
Also at z = h, 



'0, Pe=P0, 



(25) 



where po is the applied effective pressure at the 
top of the porous media, and (fio is the initial 
porosity. 

3 Non-dimensionalization 

If a length-scale d is a typical length [9] defined 
by 



d 



i^G )K = 1 + |», 
{ps ~ Pi)9 3£o 



(26) 



and the effective pressure is scaled in the fol- 
lowing way 



P 



G(p e - go) 

(Ps-PlW 



(27) 



so that p = 0(1). Here G = 1 + g is the 
value at the basin top. Compaction viscosity 
£ varies slowly with temperature as shown be- 
low in equation (32) where /3< 1 and k<1, 
and the medium viscosity rj also varies slowly 
with temperature so that the factor || does 
not change significantly because the variations 
of these two viscosities may cancel in some way 
as we now mainly focus on the region where 
temperature is relative low (< 400 K). There- 
fore, for simplicity, we take G to be constant. 



Meanwhile, we scale z with d, u s with m s , time 
t with d/rh s , permeability k with fco, and write 



T = T + |?e, 

Jo 



(28) 



where 7 is the thermal gradient, and To is the 
temperature at the basin top. we thus have 



d± d{(j)u l ) = Q 
dt dz 
dp 



(30) 



u l -u") = Xk(<f>)[^- -(!-</>)]. (31) 



dz 

The viscous relation becomes 



, du s 



p=- {1 + ( 3Q) e -^—. (32) 



where 



, _ h{ps ~ Pi) 9 _l d _ E al d 
A — : > P — ~i K — T^PT- 



pm s 



(33) 

Adding (29) and (30) together and integrating 
from the bottom, we have 



u s = -6(u l - u s ) = —u. 



(34) 



where u = (f)(u l — u s ) is the Darcy flow velocity. 
Now we have 



.dp 



(35) 



u = -\k(cf>)[-£-(l -<(>)]. (36) 

The constitutive relation for permeability k(<j)) 
is nonlinear and complicated depending on 
many parameters such as grain geometry, grain 
size distribution, materials and even the sedi- 
mentary history. For simplicity without losing 
the essence of physical mechanism of pressure 
solution concerned here, we use a simpler form 



m = (- 



(37) 



where the exponent m is derived from experi- 
mental studies. Recently, Pape et al. [17] sug- 
gested that m = 1 ~ 10 based on fractal mod- 
elling on permeability and extensive experi- 
mental studies for 640 core samples. Consid- 
ering earlier investigations [2,4,5,17], we here 
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choose a relative high value, say, m = 8, which 
is a typical value for shaly sediments. 

Different relationship of p and or u leads 
to different compaction model equations, and 
thus we have 

% = ^-<ri£-v-m, (38) 



p = A(l + /5e)e^ e £{(-f - (1 - *)]}, 



dp 



dz 4>o dz 

The boundary conditions are 
dp 



(39) 



_fl_(l_0) = O, at z = 0, (40) 



4> = <p , h = m(t)+X(- 



1-0)] at z 
(41) 

where rh(t) = O(l) is a prescribed function of 
time, which can be taken to be one for constant 
sedimentation on top of the porous media. Ob- 
viously, rh = if there is no further sedimen- 
tation and no increasing loading on top of the 
porous media. 

For simplicity, we can use a prescribed linear 
temperature profile 



e = h(t) - z. 



(42) 



It is useful for the understanding of the solu- 
tions to get an estimate for A by using val- 
ues taken from observations [6,7,9]. By using 
the typical values of pi ~ 10 3 kgm" 3 , p s ~ 
2.5 x 10 3 kg m" 3 , k ~ 10" 15 — 10" 20 m 2 , p ~ 
10" 3 Nsm 2 ,£ ~ 1 x 10 21 N s m" 2 , 
m s ~ 300 m Ma" 1 = 1 x 10 -11 m s _1 , 5 « 
10ms" 2 , G « 1, E a ~ Skcalmol" 1 , and 7 = 
0.03 Km" 1 =(30 K/1000 m); then A » 0.01 - 
-1000, p » 0.1, k » 0.2 and d » 1000 m. 

4 Numerical Simulations and 
Asymptotic Analysis 

4.1 Numerical Results 

The nonlinear diffusion equations have been 
solved by using an implicit predictor-corrector 
method. A normalized grid parameterized is 
used to get a rescaled height variable Z = 



z/h(t) in a fixed domain, which will make it 
easy to compare the results of different times 
with different values of dimensionless param- 
eters in a fixed frame. This transformation 
maps the basement of the basin to Z = and 
the basin top to Z = 1. The calculations were 
mainly implemented for the time evolutions in 
the range of t = 0.5 ~ 10 since the thickness 
in the range of 0.5km ~ 10km is the one of 
interest in the petroleum industry and in civil 
engineering. Numerical results are briefly pre- 
sented and explained below. The comparison 
with the asymptotic solutions for equilibrium 
state will be made in the next section. 

The compaction parameter A ~ 0.01 — 
— 1000, which is the ratio between the perme- 
and the sedimentation rate, defines a 
transition between the slow compaction (A << 
1) and fast compaction (A >> 1). As shown 
in [5], slow compaction is the compaction in 
a boundary layer near the basin bottom, and 
cj) (J) = 0.5, while the more interesting case 
is the fast compaction where porosity (j) reduce 
quickly. However, the effect of temperature 
gradient is not included there. Therefore, we 
now mainly investigate the effect of temper- 
ature gradient in the case of fast compaction 
when A S> 1. 

Figure 1 provides the viscous compaction 
profile of porosity versus the rescaled height 
Z = z/h(t) at different temperature gradient 
P = 0.05,0.25,0.5 for A = 100 and t = 10. 
We can see that viscous compaction profile is 
more or less parabolic in the top region. Tem- 
perature gradient greatly influence the com- 
paction behavior as pressure solution proceed, 
but the thermal effect is only of secondary im- 
portance, which is consistent with previous re- 
sults [1]. Compared with the case of constant 
permeability, porosity decreases much slower 
in the present case and this in fact implies the 
increase of the pore pressure. As the depth 
increases, the permeability k(<f>) = (0/^o) m 
may become very small as (j) < (f>o for a rel- 
ative high value of m, which will in turn con- 
strain the flow through the porous media, and 
consequently the pore fluid in sediments gets 
trapped in the lower permeability zone, result- 
ing the sudden increase of high pore pressure. 
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This can explain the general occurrance of the 
high pore pressure in sedimentary basins. 

Figure 2 gives the basin thickness h(t) as a 
function of time t for different values of A = 
0.1,10,1000. It clearly show that the moving 
boundary z = h(t) increases almost linearly 
with time t, which implies that h = const, but 
h is a function of compaction parameter A. 

To understand these phenomena and to ver- 
ify these numerical results, it would be very 
helpful if we can find some analytical solutions 
to be compared with. However, it is very diffi- 
culty to get general solutions for equations (38) 
and (39) because these equations are nonlinear 
with a moving boundary h{t). Nevertheless, it 
is still possible and very helpful to find out the 
equilibrium state and compare with the full nu- 
merical solutions [4,5]. 

4.2 Equilibrium State 

To find out the solutions for the equilibrium 
state, we must solve a nonlinear or a pair of 
nonlinear ordinary differential equations whose 
solution can usually implicitly be written in 
the quadrature form. In order to plot out and 
see the insight of the mechanism, we also need 
to solve these ordinary differential equations 
(ODEs) numerically although the solution pro- 
cedure is straightforward. However, it is prac- 
tical to get the asymptotic solutions in the ex- 
plicit form in the following cases. 

For the viscous compaction, the equilibrium 
state is governed by 

*>-<ri|;-(i-#)i> = o. 



p-*p-«»-*)i!«£ni-<i -*)]). 



where 



P = k-P. 



(43) 
(44) 



In deriving the equation (43), we have used 
the fact that j3 <C 1 and k <C 1 so that we can 
linerise the nonlinear factor in equation (32) 
by using (1 + (39) exp(-reG) « 1 - /3(h - z). 
The integration of the first equation together 



with the top boundary condition leads to 

P =^-«*-,)l|[(*ztos)]. (45) 

Subsituting this expression for p into equation 
(43) and integrating once, we obtain 



(m - h)(l - 4> Q ) 

i-<f> 



\[l-P(h-z)](- 



d 2 , i 



X {[l-f3(h-z)](m-h)(l-fo) — ( 



dz 2 y i - 



■(1-0}, 

(46) 



whose general solution can also be written in a 
quadrature. However, two distinguished limits 
are more interesting. Clearly, if A — > 0, we 
have 

h = m, 4> = 4>o, (47) 

which is the case of no compaction as discussed 
in the case of poroelastic compaction. Mean- 
while, if A — >■ oo, we have 

[l_^ ( / l _ z)]( ^_/ l)( l_ (/)o) ^ ( _l_ ) _ ( l_ 0) = o, 

(48) 

which is non-autonomous and it is difficult to 
get its general solution. However, we can as- 
sume | « 1 and perturb the above equation 
in term of /3, 



= ^°)+^W + .., 
and the leading order equation is 



(m-h)(l-4> )^{ 



8 2 , 1 



dz 2 K i - ^(o) 



which can be rewritten as 



"(1 



«.((>)> 



(49) 



= 0, 
(50) 



(m -h)(i -to)<-I = o, ^ = Y^ry 

(51) 

By using ip" = ipdip' /dip and integrating from 
h to z, we have 

(m - h)(l - <M , 2 ip J\_ 

2 ipo l-0o 

(52) 

Rearranging the above equation and changing 
variable ifi = tpo exp(^ 2 ), we get 



^0 



(i - M 



= ( Z dz. (53) 
Jh 
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After integration, we have the solution in terms 
of the original variables 0(°) and z 



i [erf 



-erf 



2 (' "*!)(*-*). 



1-0(0) i_0 o J ^(m-h) 



The first order equation is 



(54) 



<j2 jil) 

( fh _/ l )(l_^)_JP_ + ^(i) = (/j-z)^), (55) 

By using the leading order solution, the solu- 
tion for 0( x ) is simply 



«to(/i-z)-(l-0o) 



2 (1-^0) 



2(m - /i) 



x{(/i-z) 2 -2A 2 [l-™£]} (56) 
cos 



where A = y (m — — 0o). The compar- 
ison of viscous solutions (54) and (56) with 
the numerical results is shown in Figure 3 in 
the top region where the compaction profile 
is nearly at equilibrium state for A = 1000 
and t = 10 for two typical thermal gradients 
j3 = 0.1, 0.2. The agreement verifies the nu- 
merical method and the asymptotic solution 
procedure. 

4.3 Comparison With Real Data 

The numerical simulations and its comparison 
with real data are shown in Figure 4. The solid 
curve is the numerical results and real data are 
depicted by o. The real data are the borehole 
log data with a total depth of 3700 m in South 
China Sea. The rescaled height Z = z/h(t) 
varies from to 1 corresponds to a depth of 
3700 m at basement to the ocean floor. In this 
simulation, we got best fitted values of A = 
250, m = 7.3, (3 = 0.14, and t = 4.3 (or real 
time scale 14.2 Ma). 

We can see that porosity near the basin 
top decrease nearly parabolically with depth, 
and the porosity reduction only becomes 
significant at the depth h — z ~ II de- 
rived from solution (54) when its right hand 

2(1 - o )/vr(m - h)(h - z) = 0(1), that is 



n 



' 7r(m — h) 
2(l-0o) : 



(57) 



which is about 980 m for <pQ = 0.4 and m — 
h = 0.37. In other words, pressure solution 
becomes only significant at the depths greater 
than II, which is consistent with the real data. 

5 Discussions 

The present model of pressure solution in sed- 
imentary basins incorporates the effect of tem- 
perature gradient in the frame of viscous com- 
paction. Based on the pseudo-steady state 
approximations in the grain boundary diffu- 
sion process, we have been able to formulate 
a new derivation of viscous compaction rela- 
tion by using more realistic boundary condi- 
tions adjacent grain contacts. The nondimen- 
sional model equations are mainly controlled 
by two parameters A, which is the ratio of hy- 
draulic conductivity to the sedimentation rate, 
and the thermal gradient f3. Following the sim- 
ilar asymptotic analysis [5], we have been able 
to obtain the approximate solutions for either 
slow compaction (A <C 1) or fast compaction 
(A 3> 1). The more realistic and yet more in- 
teresting case is when small (but realistic) tem- 
perature gradient ^ < 1 and A > 1, and the 
equilibrium solution implies a near parabolic 
profile of porosity versus depth. Temperature 
gradient is a very large factor controlling com- 
paction process, but it is only of second impor- 
tance in the sense that it does not influence the 
parabolic shape of compaction curves since the 
shape is mainly characterized by A. However, 
for the same value of A at the same time, the 
individual curve of the compaction profile is 
essentially described by the thermal gradient. 

The numerical simulations and asymptotic 
analysis have shown that porosity-depth 
profile is near parabolic followed by a sudden 
switch of nearly uniform porosity because 
\(4>/4>o) m may become small (even A> 1 due 
to the big exponent m) at sufficiently large 
depths. In this case, the porosity profile con- 
sists of an upper part near the surface where 
the equilibrium is attained, and a lower part 
where the porosity is higher than equilibrium 
which appears to correspond accurately to nu- 
merical computations. In the near equilibrium 
region, the effect of temperature gradient is 
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very distinguished, the higher the gradient, 
the quicker the compaction proceeds. On the 
other hand, once in the nearly uniform lower 
region, the porosity is essentially uniform, the 
effect of thermal gradient is not important 
and negligible, which is consistent with pre- 
vious numerical simulations [1]. In fact, the 
permeability becomes so small that fluid gets 
trapped below this region, and compaction 
virtually stops. 
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Figure 1: Viscous compaction profile of porosity 
versus the rescaled height Z = z/h(t) at different 
temperature gradient (3 — 0.05,0.25,0.5 for A = 
100 and t = 10. The profile now is nearly parabolic. 



Figure 3: Comparison of asymptotic solutions (54) 
and (56) (dashed curves) at t = 10 for A = 1000 
with numerical results (solid curves) in the top re- 
gion (Z = z/h(t) ~ 1 or z ~ h(t)) where the 
profile is nearly at equilibrium state. The curves 
are calculated for two typical thermal gradients of 
(3 = 0.1, 0.2, and the agreement is clearly shown. 




time 



Figure 2: The basin thickness h(t) as a function 
of time t for different values of A = 0.1, 10, 1000. 
It clearly show that the moving boundary z — h(t) 
increases almost linearly with time t, which implies 
that h = const depending only on A. 




Figure 4: Comparison of numerical simulations 
(solid curves) with real borehole log data (with o). 
Z = z/h(t) is the scaled height. The best fitted 
values are A = 250, m = 7.3, /3 = 0.14, and t = 4.3. 
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